Podsumowanie

Celem poniższej analizy było znalezienie przyczyny malającej długości śledzi wyławianych w Europie. Analiza przeprowadzona została na podstawie zbioru danych z pomiarami śledzi i warunków w jakich żyli w ciągu ostatnich 60 lat. Po wyczyszczeniu danych i przeprowadzeniu analizy stwierdzono, że największy wpływ na rozmiar śledzi miała temperatura przy powierzchni wody. Im wyższa temperatura tym mniejszy śledź.

Wykorzystane biblioteki

W raporcie wykorzystano następujące biblioteki:

library(knitr)
library(DT)
library(dplyr)
library(tidyr)
library(ggplot2)
library(plotly)
library(ggcorrplot)
library(caret)

Wczytanie danych

Dane zawierają pomiary dotyczące wielkości śledzi i warunków w jakich żyją z ostatnich 60 lat.

Dane zostały zebrane z połowów komercyjnych jednostek.

df <- read.csv("./data/sledzie.csv", na.strings="?")
df <- tbl_df(df)

Opis atrybutów

Nazwa atrybutu Opis
X numer obserwacji
length długość złowionego śledzia [cm]
cfin1 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1]
cfin2 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2]
chel1 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1]
chel2 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]
lcop1 dostępność planktonu [zagęszczenie widłonogów gat. 1]
lcop2 dostępność planktonu [zagęszczenie widłonogów gat. 2]
fbar natężenie połowów w regionie [ułamek pozostawionego narybku]
recr roczny narybek [liczba śledzi]
cumf łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku]
totaln łączna liczba ryb złowionych w ramach połowu [liczba śledzi]
sst temperatura przy powierzchni wody [°C]
sal poziom zasolenia wody [Knudsen ppt]
xmonth miesiąc połowu [numer miesiąca]
nao oscylacja północnoatlantycka [mb]

Czyszczenie zbioru danych

Uzupełnienie brakujących wartości

Zbiór danych zawiera 52582 pomiarów. W zbiorze znajduje się 10094 niepełnych obserwacji, co stanowi 19.2 % wszystkich pomiarów. Liczba rekordów z wartościami NA jest zbyt duża, aby je usunąć.

Liczba brakujących wartości w poszczególnych kolumnach przedstawia się następująco:

kable(df %>% summarise_all(funs(sum(is.na(.)))))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
0 0 1581 1536 1555 1556 1653 1591 0 0 0 0 1584 0 0 0

Ze względu na to, iż wartości atrybutów często występują w podobnych grupach, to znaczy, że wartości atrybutów w sąsiadujących obserwacjach często są takie same. Brakujące wartości uzupełniono na podstawie sąsiednich rekordów - górnego lub dolnego.

na_columns <- colnames(df)[colSums(is.na(df)) > 0]
cleaned_df <- df %>% fill(na_columns, .direction="updown")

Podstawowe statystyki

Wyczyszczony zbiór danych składa się z 52582 wierszy (obserwacji) i 16 kolumn (atrybutów).

kable(summary(cleaned_df %>% select(X:lcop2)))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2
Min. : 0 Min. :19.0 Min. : 0.0000 Min. : 0.0000 Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849
1st Qu.:13145 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778 1st Qu.: 2.469 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808
Median :26291 Median :25.5 Median : 0.1111 Median : 0.7012 Median : 5.750 Median :21.435 Median : 7.0000 Median :24.859
Mean :26291 Mean :25.3 Mean : 0.4457 Mean : 2.0255 Mean :10.003 Mean :21.215 Mean : 12.8079 Mean :28.419
3rd Qu.:39436 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936 3rd Qu.:11.500 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232
Max. :52581 Max. :32.5 Max. :37.6667 Max. :19.3958 Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736 
kable(summary(cleaned_df %>% select(fbar:nao)))
fbar recr cumf totaln sst sal xmonth nao
Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137 Min. :12.77 Min. :35.40 Min. : 1.000 Min. :-4.89000
1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068 1st Qu.:13.60 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
Median :0.3320 Median : 421391 Median :0.23191 Median : 539558 Median :13.86 Median :35.51 Median : 8.000 Median : 0.20000
Mean :0.3304 Mean : 520367 Mean :0.22981 Mean : 514973 Mean :13.88 Mean :35.51 Mean : 7.258 Mean :-0.09236
3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351 3rd Qu.:14.16 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595 Max. :14.73 Max. :35.61 Max. :12.000 Max. : 5.08000

Analiza danych

Rozkład wartości atrybutów

p <- ggplot(gather(cleaned_df %>% select(-X, -xmonth)), aes(x = value)) + 
  geom_histogram() + 
  facet_wrap(~key, scales="free") + 
  theme_light()

ggplotly(p)

Zmiana poszczególnych atrybutów w czasie

Poniżej zaprezentowano wykresy zmiany poszczególnych atrybutów na przestrzeni kolejnych obserwacji. Zmianę atrybutów w czasie zaprezentowano za pomocą linii trendu.

plankton_df <- cleaned_df %>%
  select(X, cfin1:lcop2) %>%
  gather(plankton_name, plankton_density, cfin1:lcop2)

p <- ggplot(plankton_df, aes(x=X, y=plankton_density, color = plankton_name)) + 
  geom_smooth() + 
  labs(title = "Zmiana dostępności planktonu", x="Nr obserwacji", 
       y="Zagęszczenie planktonu") + 
  scale_color_discrete(name=NULL) +
  theme_light()

ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=fbar)) +
  geom_smooth() +
  labs(title = "Zmiana natężenia połowów w regionie", x="Nr obserwacji", 
       y="Natężenie połowów w regionie [ułamek pozostawionego narybku]") + 
  theme_light()

ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=recr)) +
  geom_smooth() +
  labs(title = "Zmiana rocznego narybku", x="Nr obserwacji", 
       y="Roczny narybek [liczba śledzi]") + 
  scale_y_continuous(labels=scales::comma) + 
  theme_light()

ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=cumf)) +
  geom_smooth() +
  labs(title = "Zmiana łącznego rocznego natężenia połowów w regionie", 
       x="Nr obserwacji", y="Ułamek pozostawionego narybku") + 
  theme_light()

ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=totaln)) +
  geom_smooth() +
  labs(title = "Zmiana łącznej liczby ryb złowionych w ramach połowu", 
       x="Nr obserwacji", y="Liczba śledzi") + 
  scale_y_continuous(labels=scales::comma) + 
  theme_light()

ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=sst)) +
  geom_smooth() +
  labs(title = "Zmiana temperatury przy powierzchni wody", x="Nr obserwacji", 
       y="Temperatura [°C]") + 
  theme_light()

ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=sal)) +
  geom_smooth() +
  labs(title = "Zmiana poziomu zasolenia wody", x="Nr obserwacji", 
       y="Zasolenie [Knudsen ppt]") + 
  theme_light()

ggplotly(p)
p <- ggplot(cleaned_df, aes(x=X, y=nao)) +
  geom_smooth() +
  labs(title = "Zmiana oscylacji północnoatlantyckiej", x="Nr obserwacji", 
       y="Oscylacja północnoatlantycka [mb]") + 
  theme_light()

ggplotly(p)

Zmiana rozmiaru śledzia w czasie

Zmianę wielkości śledzia zaprezentowano na interaktywnym wykresie za pomocą linii trendu. Oś x oznacza liczbę porządkową - numer obserwacji, natomiast oś y rozmiar śledzia w cm.

p <- ggplot(cleaned_df, aes(x=X, y=length)) + 
  geom_smooth() + 
  labs(x="Nr obserwacji", y="Długość [cm]") + 
  theme_light()

ggplotly(p)

Korelacja pomiędzy atrybutami

Korelacje pomiędzy poszczególnymi atrybutami przedstawiono na interaktywnej mapie korelacji. Do obliczenia współczynnika korelacji wykorzystano metodę Pearsona.

Aby wyświetlić nazwy skorelowanych zmiennych i wartości korelacji należy najechać kursorem na poszczególne komórki na wykresie.

cor_matrix <- cleaned_df %>% 
  select(-X) %>% 
  cor(use = "all.obs", method="pearson")

corr_plot <- cor_matrix %>% 
  ggcorrplot(type="lower", legend.title = "Współczynnik Pearsona") +
  labs(x = 'Atrybut 1', y = 'Atrybut 2') + 
  theme_light() +
  theme(axis.ticks = element_blank()) 

ggplotly(corr_plot)

Na podstawie wykresu prezentującego korelację pomiędzy atrybutami można zauważyć, że atrybuty lcop1 i chel1 są ze sobą najsilniej skorelowane. Atrybut length jest najsiliniej dodatnio skorelowany z atrybutami fbar oraz lcop1, natomiast najsilniej skorelowany ujemnie z atrybutami sst oraz nao.

Regresor przewidujący rozmiar śledzia

Redukcja korelacji

W celu redukcji korelacji między atrybutami wykorzystano funkcję ‘findCorrelation’ z pakietu ‘caret’ z parametrem cutoff równym 0.8, która zwraca atrybuty do usunięcia.

attributes_to_remove <- cor_matrix %>% findCorrelation(cutoff = 0.8, names = TRUE)

Atrybuty wybrane do usunięcia: lcop2, cumf, chel1.

Trenowanie modelu

Do predykcji pominięto atrybuty lcop2, cumf, chel1 oraz X. Zbiór danych podzielono na zbiór uczący i testowy w proporcjach 75/25. Zbiór walidacyjny został utworzony przy użyciu wielokrotnej oceny krzyżowej (ang. repeated cross-validation) z liczbą podziałów równą 2 i liczbą powtórzeń równą 5.

in_training_data <- createDataPartition(y = cleaned_df$length, p = 0.75, list = FALSE)

training_data <- cleaned_df[in_training_data, ] %>% select(-c(X, attributes_to_remove))
testing_data <- cleaned_df[-in_training_data, ]

ctrl <- trainControl(method = "repeatedcv", number = 2, repeats = 5)

Poniższy wykres przedstawia podobieństwo rozkładów danych treningowych i testowych.

ggplot() +
  geom_density(aes(length, fill = "Treningowy"), training_data, alpha = 0.6) +
  geom_density(aes(length, fill = "Testowy"), testing_data, alpha = 0.6) +
  labs(x = "Długość [cm]", y = "Gęstość", fill = "Zbiór danych") +
  theme_light()

Do przewidywania długości śledzia wykorzystano metodę Random Forest. Do oceny dokładności predykcji wykorzystano miary R2 i RMSE.

fit <- train(length ~ .,
             data = training_data,
             method = "rf",
             trControl = ctrl,
             ntree = 10)

fit
## Random Forest 
## 
## 39438 samples
##    11 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 5 times) 
## Summary of sample sizes: 19720, 19718, 19719, 19719, 19719, 19719, ... 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE      
##    2    1.166928  0.5021753  0.9225294
##    6    1.148083  0.5181004  0.9045890
##   11    1.151060  0.5160672  0.9054773
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 6.

Wartości miar dla zbioru testowego:

prediction <- predict(fit, testing_data)

post_resample <- postResample(pred = prediction,
                              obs = testing_data$length)
post_resample
##      RMSE  Rsquared       MAE 
## 1.1372567 0.5251405 0.8957069

Poniższy wykres przedstawia wartości zbioru testowego oraz wartości przewidziane przez regresor.

prediction_comparison_df <- tibble(X = testing_data$X, 
                                   actual = testing_data$length,
                                   predicted = prediction)

ggplot(prediction_comparison_df, aes(x = X)) +
  geom_smooth(aes(y = actual, color = "Wartość rzeczywista")) +
  geom_smooth(aes(y = predicted, color = "Wartość przewidziana")) +
  labs(color = "Wartości", x = "Nr obserwacji", y = "Długość [cm]") +
  theme_light()

Ważność atrybutów

ggplot(varImp(fit)) + 
  labs(x = "Atrybut", y = "Ważność") +
  theme_light()

Patrząc na powyższy wykres okazało się, że w przewidywaniu rozmiaru śledzia najważniejsza była temperatura przy powierzchni wody (sst). Mniejszy wpływ miały również natężenie połowów w regionie (fbar) oraz roczny narybek (recr).

Na podstawie całej analizy można stwierdzić, że na rozmiar śledzia największy wpływ miała temperatura przy powierzchni wody. Im wyższa temperatura tym mniejszy śledź.